import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
import datetime as dt
import folium
import random
from folium.plugins import MeasureControl
from folium import plugins
from pulp import *
import json
pd.set_option("display.max_columns", 100)
pd.set_option("display.max_rows",700)
from haversine import haversine

data = pd.read_csv("Boston Crimes Data 2017.csv")

data.head(10)

	INCIDENT_NUMBER	OFFENSE_CODE	OFFENSE_CODE_GROUP	OFFENSE_DESCRIPTION	DISTRICT	REPORTING_AREA	SHOOTING	OCCURRED_ON_DATE	YEAR	MONTH	DAY_OF_WEEK	HOUR	UCR_PART	STREET	Lat	Long	Location
0	I172067337	3410	Towed	TOWED MOTOR VEHICLE	C6	202	NaN	2017-08-16 01:35:00	2017	8	Wednesday	1	Part Three	W SECOND ST	42.340643	-71.050031	(42.34064258, -71.05003074)
1	I172075354	3006	Medical Assistance	SICK/INJURED/MEDICAL - PERSON	B2	319	NaN	2017-09-11 10:13:00	2017	9	Monday	10	Part Three	SEAVER ST	42.306040	-71.086212	(42.30603974, -71.08621240)
2	I172050693	3114	Investigate Property	INVESTIGATE PROPERTY	A1	112	NaN	2017-06-24 13:28:00	2017	6	Saturday	13	Part Three	AVENUE DE LAFAYETTE	42.353320	-71.061023	(42.35331987, -71.06102314)
3	I172019883	3115	Investigate Person	INVESTIGATE PERSON	C11	342	NaN	2017-03-13 16:15:00	2017	3	Monday	16	Part Three	HOMES AVE	42.302884	-71.069188	(42.30288397, -71.06918798)
4	I172104661	3114	Investigate Property	INVESTIGATE PROPERTY	A1	81	NaN	2017-12-07 00:00:00	2017	12	Thursday	0	Part Three	COMMERCIAL ST	42.368305	-71.056188	(42.36830502, -71.05618848)
5	I172018265	3125	Warrant Arrests	WARRANT ARREST	D4	165	NaN	2017-03-07 18:25:00	2017	3	Tuesday	18	Part Three	HARRISON AVE	42.342850	-71.065162	(42.34285014, -71.06516235)
6	I172016666	3018	Medical Assistance	SICK/INJURED/MEDICAL - POLICE	B2	182	NaN	2017-03-01 20:24:00	2017	3	Wednesday	20	Part Three	ALBION ST	42.322038	-71.071404	(42.32203840, -71.07140396)
7	I172023186	2610	Other	TRESPASSING	D4	273	NaN	2017-03-25 04:44:00	2017	3	Saturday	4	Part Two	WASHINGTON ST	42.335435	-71.078402	(42.33543510, -71.07840204)
8	I172000228	802	Simple Assault	ASSAULT SIMPLE - BATTERY	C11	364	NaN	2017-01-01 14:00:00	2017	1	Sunday	14	Part Two	LESLIE ST	42.295147	-71.058608	(42.29514664, -71.05860832)
9	I172068767	617	Larceny	LARCENY THEFT FROM BUILDING	A1	77	NaN	2017-08-20 14:27:00	2017	8	Sunday	14	Part One	CANAL ST	42.364895	-71.060502	(42.36489527, -71.06050249)

most_freq_crimes = dict(data.groupby("OFFENSE_CODE_GROUP")["INCIDENT_NUMBER"].count().sort_values(ascending=False))
print(dict(list(most_freq_crimes.items())[:5]))

{'Motor Vehicle Accident Response': 6022, 'Larceny': 3922, 'Medical Assistance': 3862, 'Investigate Person': 3252, 'Other': 2657}

motor_vehicle_accidents_df = data[data['OFFENSE_CODE_GROUP']=='Motor Vehicle Accident Response'].reset_index(drop=True)
location_data = motor_vehicle_accidents_df[["INCIDENT_NUMBER","DISTRICT","YEAR","Lat","Long","Location"]]
location_data=location_data.dropna(subset = ['Lat','DISTRICT'])
location_data.isnull().sum()

INCIDENT_NUMBER    0
DISTRICT           0
YEAR               0
Lat                0
Long               0
Location           0
dtype: int64

location_data["Location"]=location_data["Location"].apply(lambda x:x.strip("()").strip(" "))
location_data["Location"]=location_data["Location"].apply(lambda x:str(x.split(" ")[0]+x.split(" ")[1]))
location_data=location_data[~location_data.Long.isin([-1])].reset_index(drop=True)
location_data=location_data[location_data.DISTRICT.isin(['E13','B3','C11','B2','D4','C6','A1'])].reset_index(drop=True)

my_colours = {'C11': 'purple',
 'B3': 'orange',
 'E13': 'red',
 'A1': 'yellow',
 'C6': 'green',
 'B2': 'blue',
 'D4': 'pink'}


location_data["my_colours"] = location_data["DISTRICT"].apply(lambda x:my_colours[x])
location_data = location_data.drop_duplicates(subset="INCIDENT_NUMBER", keep="first").reset_index(drop=True)

a = location_data.sample(2000).reset_index(drop=True)
WHS_COORD = [location_data["Lat"][0],location_data["Long"][0]]
map_nyc = folium.Map(location=WHS_COORD, zoom_start=12, width=740, height=500)

for i in range(len(a)):
    folium.CircleMarker([a["Lat"][i],a["Long"][i]], radius=5,color=a["my_colours"][i],
                        popup=a["DISTRICT"][i]).add_to(map_nyc)
    
plugins.Fullscreen(position='topleft').add_to(map_nyc)
map_nyc.add_child(MeasureControl())

map_nyc

Make this Notebook Trusted to load map: File -> Trust Notebook

optimization_data = location_data[["INCIDENT_NUMBER","Lat","Long",'my_colours']]

Potential_Station_Number=300
Station_ids=["ST"+str(x) for x in range(1,Potential_Station_Number+1)]
randomlist = random.sample(range(0, optimization_data.index.max()), len(Station_ids))
potential_stations=pd.DataFrame({"Station_id":Station_ids,
              "Station_Lat":optimization_data.loc[randomlist,"Lat"].values.tolist(),
             "Station_Long":optimization_data.loc[randomlist,"Long"].values.tolist()})

WHS_COORD = [optimization_data["Lat"][0],optimization_data["Long"][0]]
map_nyc = folium.Map(location=WHS_COORD, zoom_start=12, width=740, height=500)


for i in range(len(optimization_data)):
    
    folium.CircleMarker([optimization_data["Lat"][i],optimization_data["Long"][i]], 
                        radius=5,color=optimization_data["my_colours"][i], fill_color='#0080bb',
                        popup={optimization_data["INCIDENT_NUMBER"][i]}).add_to(map_nyc)
    

for j in range(len(potential_stations)):
    folium.Marker([potential_stations["Station_Lat"][j],potential_stations["Station_Long"][j]], radius=5,
                  popup=potential_stations["Station_id"][j]).add_to(map_nyc) 
map_nyc

Make this Notebook Trusted to load map: File -> Trust Notebook

Data Preperation for Optimization Algorithm

N = optimization_data["INCIDENT_NUMBER"].nunique()

print("The dimension : " + str(N)) 
res = [list(range(1 + N * i, 1 + N * (i + 1))) 
                            for i in range(N)]  
res=pd.DataFrame(res)

res=res[:len(potential_stations)]
res.columns=optimization_data["INCIDENT_NUMBER"]
res[res > 0] = 0
res.index=potential_stations["Station_id"].values

potential_coords=[]
for j in range(len(potential_stations)):
    potential_coords.append(str(potential_stations["Station_Lat"][j])+","+str(potential_stations["Station_Long"][j]))

demand_coords=[]
for j in range(len(optimization_data)):
    demand_coords.append(str(optimization_data["Lat"][j])+","+str(optimization_data["Long"][j]))


res["potential_coords"] = potential_coords
demand_coords.append("0")
demand_coords=pd.DataFrame(demand_coords).T
demand_coords.columns=res.columns
res=res.append(demand_coords)

last_data=pd.DataFrame()
my_list=[]
for j in range(0,(len(res.columns)-1)):
    for i in range(len(res)-1):
        my_length=int(haversine([float(res.loc[res.index[i]]["potential_coords"].split(",")[0]),
              float(res.loc[res.index[i]]["potential_coords"].split(",")[1])],
             [float(res[res.columns[j]][0].split(",")[0]),float(res[res.columns[j]][0].split(",")[1])])*1000)
        my_list.append(my_length)  
    last_data[res.columns[j]] = my_list
    my_list=[]
last_data.index=potential_stations["Station_id"].values

The dimension : 3528

/var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/ipykernel_47073/350201403.py:36: PerformanceWarning: DataFrame is highly fragmented.  This is usually the result of calling `frame.insert` many times, which has poor performance.  Consider joining all columns at once using pd.concat(axis=1) instead.  To get a de-fragmented frame, use `newframe = frame.copy()`
  last_data[res.columns[j]] = my_list

last_data.head()

	I172034460	I172026267	I172002734	I172055514	I172025030	I172003693	I172060524	I172015101	I172035682	I172028199	I172054480	I172065426	I172039140	I172018430	I172021387	I172078710	I172012419	I172062371	I172066844	I172092769	I172017788	I172012405	I172080591	I172068926	I172010516	I172091713	I172104960	I172072011	I172028902	I172078242	I172021745	I172095142	I172020528	I172046176	I172047611	I172095279	I172089232	I172083138	I172039512	I172083893	I172008642	I172016632	I172025006	I172044610	I172036642	I172082738	I172075413	I172036206	I172009578	I172021561	...	I172027418	I172107075	I172073840	I172047408	I172084997	I172087115	I172065757	I172098580	I172015974	I172107439	I172053387	I172106757	I172008873	I172084857	I172092971	I172015767	I172010184	I172052079	I182000005	I172103242	I172070059	I172019465	I172009946	I172073888	I172102986	I172012979	I172018466	I172039985	I172052449	I172002127	I172087522	I172008116	I172026815	I172062951	I172063987	I172035554	I172015446	I172081504	I172066094	I172011629	I172106209	I172044243	I172027109	I172024106	I172002830	I172022354	I172007363	I172062252	I172070810	I172106706
ST1	3178	1622	5277	1761	1836	680	1725	6044	3935	2007	5902	3873	4147	6767	2837	2032	2213	3305	3730	7081	6648	6752	4052	5463	3227	5630	6172	4350	7841	1152	2189	3934	5860	7129	2820	4162	5898	1479	4920	2513	1701	1932	6294	4371	2408	3262	2382	7263	5892	7101	...	5569	822	8324	6061	1154	7893	1250	5196	4408	1154	2456	5953	2691	3811	8826	3375	6419	8318	1918	1503	2545	5098	2199	4422	7035	4792	2212	719	6996	7081	1865	1049	4874	2114	3678	2095	5783	5530	1741	1748	2690	3068	3832	5656	3178	1822	7483	1700	2795	2545
ST2	812	5287	4513	5468	2034	3257	2968	2176	3026	2226	2398	3849	4712	4489	4255	5767	2296	5735	1976	3435	3532	3687	4735	2265	2891	2962	2531	1445	4428	3751	3945	3436	4294	4226	1173	1607	2168	2395	4162	2031	2199	3633	4305	2556	6146	654	2074	3733	3666	3397	...	1722	4469	4955	4552	3807	4206	2634	2290	4553	3807	1416	3490	1804	3408	5096	2707	2654	4923	3705	3196	3509	2644	4470	4536	3778	3999	6083	4120	4547	3435	2292	3445	1097	1856	426	2367	1912	1702	2131	5367	1889	2157	74	4705	812	2052	5633	2186	4079	2970
ST3	3668	1898	6833	3324	2445	1669	3286	6114	5322	1854	6420	5506	5932	7923	4682	3505	3454	5222	4784	7493	7421	7556	5853	6147	4675	6571	5770	5070	8452	788	335	5445	7219	8027	2893	4995	6226	2318	6453	3641	2559	439	7535	5517	3816	3412	3517	7780	7026	6689	...	5497	2648	8971	7449	776	8262	2268	6004	6125	776	2914	6998	3713	5334	9135	4742	6701	8946	363	3164	959	6088	426	6133	7738	6306	3143	1301	8092	7493	1827	2813	5156	2305	3705	3361	5830	5421	2436	1878	2236	4261	3994	7184	3668	2513	8818	2229	4612	4085
ST4	3516	1186	6042	2172	2142	945	2457	6271	4640	1971	6307	4653	4960	7414	3655	2380	2831	4052	4311	7456	7157	7273	4867	5927	3945	6193	6205	4813	8297	543	1491	4679	6569	7687	2961	4666	6228	1858	5677	3094	2115	1318	6967	4996	2718	3459	2962	7680	6526	7140	...	5732	1465	8796	6780	509	8256	1693	5705	5210	509	2744	6556	3235	4558	9172	4067	6737	8784	1281	2266	1972	5674	1401	5222	7521	5543	2209	115	7624	7456	1862	1839	5171	2252	3839	2718	5998	5677	2079	1261	2594	3690	4055	6414	3516	2167	8190	1951	3610	3285
ST5	4157	1831	7284	3448	2933	2092	3717	6566	5789	2340	6904	5941	6332	8410	5057	3594	3922	5486	5270	7968	7911	8047	6248	6636	5133	7062	6166	5561	8936	1183	625	5899	7693	8518	3373	5486	6699	2795	6907	4118	3039	900	8018	6001	3873	3891	3993	8260	7513	7075	...	5939	2890	9456	7921	1153	8732	2731	6495	6545	1153	3405	7488	4197	5786	9599	5209	7169	9431	829	3577	1249	6578	65	6554	8226	6762	3104	1554	8582	7968	2317	3199	5630	2795	4172	3826	6282	5859	2921	1774	2689	4739	4468	7640	4157	2999	9296	2717	4998	4528

5 rows × 3528 columns

OPTIMIZATION MODELLING AND SOLVING

my_range = range(3,13,3)
my_analysis_result = {}
for my_st_number in my_range:
    potentials=potential_stations["Station_id"].values.tolist()#[station:station+1]
    demands=optimization_data["INCIDENT_NUMBER"].values.tolist()
    distance=last_data.to_dict('index')
    prob = LpProblem("Transportation", LpMinimize)
    routes  =[(i,j) for i in potentials for j in demands]
    
    #DECISION VARIABLES
    amount_vars = LpVariable.dicts("X",(potentials,demands),lowBound=0, upBound=1, cat='Binary')
    station_vars = LpVariable.dicts("St",potentials,lowBound=0, upBound=1, cat='Binary')
    
    # OBJECTIVE
    
    prob += lpSum(amount_vars[i][j]*distance[i][j] for (i,j) in routes)

    # CONSTRAINTS
    
    for j in demands:
        prob += lpSum(amount_vars[i][j] for i in potentials) == 1

    for i in demands:
        for j in potentials:
            prob += amount_vars[j][i] <= 1*station_vars[j]

    prob += lpSum(station_vars[i] for i in potentials) == my_st_number

    
    # SOLVE THE PROBLEM
    prob.solve()
    
    
    
    wh_loc=[]
    hex_loc=[]
    assign=[]
    for v in prob.variables()[len(potentials):]:
        wh_loc.append(v.name.split("_")[1])
        hex_loc.append(v.name.split("_")[2])
        assign.append(v.varValue)

    df=pd.DataFrame({"Station_id":wh_loc,"hex_loc":hex_loc,"assign":assign})
    df_ones = df[df["assign"]!=0].reset_index(drop=True)
    my_data = df_ones
    my_data_whs = my_data.drop_duplicates("Station_id").reset_index(drop=True)
    my_data_whs = my_data_whs.merge(potential_stations,left_on="Station_id",right_on="Station_id",how="left")
    my_data_whs = my_data_whs.drop(columns=['hex_loc'])
    
    my_analysis_result["{} stations".format(my_st_number)] = my_data_whs.to_dict()

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/sdemir/opt/anaconda3/lib/python3.9/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/80ed06c560d447209728efab36733782-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/80ed06c560d447209728efab36733782-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 1061934 COLUMNS
At line 7412557 RHS
At line 8474487 BOUNDS
At line 9533188 ENDATA
Problem MODEL has 1061929 rows, 1058700 columns and 3175500 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 5.47354e+06 - 3339.62 seconds
Cgl0004I processed model has 1061929 rows, 1058700 columns (1058700 integer (1058700 of which binary)) and 3175500 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 5.47354e+06
Cbc0038I Before mini branch and bound, 1058700 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (3529.41 seconds)
Cbc0038I After 3529.44 seconds - Feasibility pump exiting with objective of 5.47354e+06 - took 1.88 seconds
Cbc0012I Integer solution of 5473545 found by feasibility pump after 0 iterations and 0 nodes (3529.80 seconds)
Cbc0001I Search completed - best objective 5473545, took 0 iterations and 0 nodes (3530.98 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 5.47354e+06 to 5.47354e+06
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                5473545.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             3514.39
Time (Wallclock seconds):       3534.35

Option for printingOptions changed from normal to all
Total time (CPU seconds):       3520.23   (Wallclock seconds):       3540.68

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/sdemir/opt/anaconda3/lib/python3.9/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/78e4448de54e4581aa7ee78a50cf892a-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/78e4448de54e4581aa7ee78a50cf892a-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 1061934 COLUMNS
At line 7412557 RHS
At line 8474487 BOUNDS
At line 9533188 ENDATA
Problem MODEL has 1061929 rows, 1058700 columns and 3175500 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 3.90138e+06 - 2209.22 seconds
Cgl0004I processed model has 1061929 rows, 1058700 columns (1058700 integer (1058700 of which binary)) and 3175500 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 3.90138e+06
Cbc0038I Before mini branch and bound, 1058700 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (2390.25 seconds)
Cbc0038I After 2390.29 seconds - Feasibility pump exiting with objective of 3.90138e+06 - took 1.86 seconds
Cbc0012I Integer solution of 3901380 found by feasibility pump after 0 iterations and 0 nodes (2390.63 seconds)
Cbc0001I Search completed - best objective 3901380, took 0 iterations and 0 nodes (2391.86 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 3.90138e+06 to 3.90138e+06
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                3901380.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             2384.48
Time (Wallclock seconds):       2395.42

Option for printingOptions changed from normal to all
Total time (CPU seconds):       2390.34   (Wallclock seconds):       2401.76

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/sdemir/opt/anaconda3/lib/python3.9/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/04edee90ff9e45e8a7c81e599d2b4447-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/04edee90ff9e45e8a7c81e599d2b4447-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 1061934 COLUMNS
At line 7412557 RHS
At line 8474487 BOUNDS
At line 9533188 ENDATA
Problem MODEL has 1061929 rows, 1058700 columns and 3175500 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 3.10466e+06 - 2173.89 seconds
Cgl0004I processed model has 1061929 rows, 1058700 columns (1058700 integer (1058700 of which binary)) and 3175500 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 3.10466e+06
Cbc0038I Before mini branch and bound, 1058700 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (2353.44 seconds)
Cbc0038I After 2353.48 seconds - Feasibility pump exiting with objective of 3.10466e+06 - took 1.73 seconds
Cbc0012I Integer solution of 3104665 found by feasibility pump after 0 iterations and 0 nodes (2353.88 seconds)
Cbc0001I Search completed - best objective 3104665, took 0 iterations and 0 nodes (2354.97 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 3.10466e+06 to 3.10466e+06
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                3104665.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             2347.74
Time (Wallclock seconds):       2358.42

Option for printingOptions changed from normal to all
Total time (CPU seconds):       2353.59   (Wallclock seconds):       2364.74

Welcome to the CBC MILP Solver 
Version: 2.10.3 
Build Date: Dec 15 2019 

command line - /Users/sdemir/opt/anaconda3/lib/python3.9/site-packages/pulp/apis/../solverdir/cbc/osx/64/cbc /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/f4f8273d2af84135836838467443c068-pulp.mps timeMode elapsed branch printingOptions all solution /var/folders/v4/sgw2lcxd7lq92jq9np2lwppm0000gq/T/f4f8273d2af84135836838467443c068-pulp.sol (default strategy 1)
At line 2 NAME          MODEL
At line 3 ROWS
At line 1061934 COLUMNS
At line 7412557 RHS
At line 8474487 BOUNDS
At line 9533188 ENDATA
Problem MODEL has 1061929 rows, 1058700 columns and 3175500 elements
Coin0008I MODEL read with 0 errors
Option for timeMode changed from cpu to elapsed
Continuous objective value is 2.62176e+06 - 185.53 seconds
Cgl0004I processed model has 1061929 rows, 1058700 columns (1058700 integer (1058700 of which binary)) and 3175500 elements
Cutoff increment increased from 1e-05 to 0.9999
Cbc0038I Initial state - 0 integers unsatisfied sum - 0
Cbc0038I Solution found of 2.62176e+06
Cbc0038I Before mini branch and bound, 1058700 integers at bound fixed and 0 continuous
Cbc0038I Mini branch and bound did not improve solution (360.03 seconds)
Cbc0038I After 360.06 seconds - Feasibility pump exiting with objective of 2.62176e+06 - took 1.70 seconds
Cbc0012I Integer solution of 2621763 found by feasibility pump after 0 iterations and 0 nodes (360.42 seconds)
Cbc0001I Search completed - best objective 2621763, took 0 iterations and 0 nodes (361.52 seconds)
Cbc0035I Maximum depth 0, 0 variables fixed on reduced cost
Cuts at root node changed objective from 2.62176e+06 to 2.62176e+06
Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Gomory was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Knapsack was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Clique was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
MixedIntegerRounding2 was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
FlowCover was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
TwoMirCuts was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
ZeroHalf was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)

Result - Optimal solution found

Objective value:                2621763.00000000
Enumerated nodes:               0
Total iterations:               0
Time (CPU seconds):             358.55
Time (Wallclock seconds):       364.70

Option for printingOptions changed from normal to all
Total time (CPU seconds):       364.49   (Wallclock seconds):       371.13

optimization_data.to_csv('optimization/original_boston_data.csv',index=False)

import pickle
my_analysis_result['potential_stations'] = potential_stations.to_dict()

with open('optimization/optimal_stations.pickle', 'wb') as handle:
    pickle.dump(my_analysis_result, handle, protocol=pickle.HIGHEST_PROTOCOL)